Atelier Statistiques Avancées
28e Colloque du Chapitre Saint-Laurent
1. Introduction
Contexte
Ce tutoriel a été développé dans le cadre du 28e colloque du Chapitre Saint-Laurent, présenté le 29 Mai 2024 à l’Université Laval, Québec, Canada.
Nous verrons les concepts de base qui vous permettront d’implémenter des modèles hiérarchiques avec le logiciel R, et d’interpréter les résultats en utilisant l’inférence Bayésienne.
Les objectifs d’apprentissage sont:
- De se familiariser avec la théorie sur les modèles linéaires hiérarchiques
- De reconnaître les variantes des modèles linéaires ((H)(G)LM)
- D’implémenter un modèle hiéarchique avec R
- D’interpréter les résultats avec l’inférence Bayésienne
librairies
Voici la liste des librairies que nous allons utiliser pour l’atelier. Assurez-vous de les avoir téléchargées.
# Pour manipuler les données
library(data.table)
# Pour estimer les paramètres des modèles
library(brms)
# Pour diagnostiquer les modèles bayésiens
library(bayesplot)
# Pour nos graphiques
library(ggplot2)
library(ggpubr)
library(viridis)
# Pour générer de jolies tables
library(knitr)Jeu de données
Le jeu de données que nous utiliserons pour l’atelier se nomme ADORE et provient de l’article suivant:
Vous pouvez trouver les données originales sur le portail du Eawag Research Data Institutional Collection.
Le jeu de données ADORE contient une myriade d’informations sur des expériences écotoxicologiques effectuées sur des espèces de poisson, de crustacés, et d’algues. L’objectif des données ADORE est de concerter des informations sur la toxicité aquatique accrue de différents composés chimiques sur des espèces pour faire de la modélisation prédictive.
2. Une note sur l’inférence Bayésienne
L’inférence Bayésienne est une toute autre façon de penser pour la réalisation et l’interprétation d’analyses statistiques. Cela vous apparaîtra comme assez différent de ce à quoi vous êtes habitué.e.s si vous avez appris les statistiques classiques (fréquentistes).
Le théorème de Bayes
En inférence Bayésienne, on se sert de nos connaissances à priori ainsi que des données pour faire de l’inférence. Nous pouvons représenter ceci avec le théorème de Bayes qui va comme suit:
\[P(H \mid E) = \frac{P(H) \ P(E \mid H)}{P(E)}\]
Pour résumer l’équation:
- \(P(H)\) est la probabilité d’observer notre hypothèse, c’est la distribution à priori
- \(P(E)\) est la probabilité d’observer de nouvelles preuves sous toutes les hypothèses possibles
- \(P(H \mid E)\) est la probabilité d’observer notre hypothèse selon l’évidence, c’est la distribution à postériori
- \(P(E \mid H)\) est la probabilité d’observer l’évidence si notre hypothèse est vraie, c’est la distribution de vraisemblance
Plus intuitivement, on peut aussi l’illustrer comme suit:
Illustration du théorème de Bayes issue de Yanagisawa et al. 2019
Ainsi, l’inférence Bayésienne est une approche où l’on adapte nos connaissances (à posteriori) à partir de nos hypothèses (connaissances à priori) et les preuves que nous avons.
La distribution à posteriori est donc l’élément central que nous cherchons à estimer. Rappelez-vous que chaque fois que vous voyez le terme “à posteriori” dans ce tutoriel, nous faisons référence à la distribution d’un paramètre d’intérêt estimé par le modèle.
En prime, voici un lien vers une application qui vous permet de jouer avec ces trois distributions pour mieux comprendre comment les connaissances à priori et les preuves modifient la distribution à posteriori:
Méthode d’estimation
Il existe de nombreuses façons d’approximer la distribution à posteriori. Dans ce tutoriel, nous utiliserons la méthode de Monte Carlo Hamiltonienne (HMC) et l’échantillonneur No-U-Turn (NUTS) pour estimer la distribution à posteriori. Ce sont des algorithmes spécifiques dans la grande famille d’estimation de Monte Carlo par chaîne de Markov (MCMC).
Pour ce faire nous effectuerons tous nos modèles avec la librairie brms, une interface R pour le langage probabiliste STAN. La librairie brms utilise l’algorithme NUTS pour estimer les paramètres. Je vous recommande vivement de lire les vignettes de la librairie pour mieux comprendre les capacités et le fonctionnement de brms.
Les avantages
L’inférence Bayésienne comporte certains avantages par rapport à l’approche statistique fréquentiste.
- Elle met l’accent sur la distribution et l’incertitude au lieu d’estimés uniques
- Elle permet une interprétation plus naturelle (probabiliste)
- On peut facilement critiquer les modèles
Pour aller plus loin
La théorie sous-jacente au théorème de Bayes dépasse largement le cadre de ce tutoriel, donc veuillez consulter les ressources pertinentes dans la section Références. Une très bonne introduction que je ne peux recommander assez est l’excellent livre “Statistical Rethinking” de Richard McElreath. Richard McElreath publie également l’ensemble du contenu du livre sous forme de vidéos sur YouTube.
3. Structure des données
Description générale
Pour l’atelier, nous utiliserons une version des données nettoyées et prêtes pour l’analyse qui se trouvent sur ce dépôt GitHub dans le dossier data. Ce dépôt a été préparé spécifiquement pour l’atelier.
Voici les spécifications:
- Le jeu de données servira à faire des modèles linéaires généralisés (GLM) et des modèles linéaires généralisés hiérarchiques (HGLM) univariés
- Nous appliquerons les modèles aux données d’animaux seulement. Nous appelerons le jeu de données
anim
Importons les données à partir du dépôt GitHub de l’atelier:
github <- "https://raw.githubusercontent.com/quantitative-ecologist"
depot <- "Atelier-CSL-EcotoQ"
dossier <- "main/data"
anim <- fread(
file.path(github, depot, dossier, "animal_dat.csv"),
stringsAsFactors = TRUE
)Variables
Une description détaillée des variables du jeu de données est présentée en Annexe de l’article de Schür et al. (2023).
Dans le cadre de cet atelier, nous utiliserons un échantillon de ces variables.
Variables sur les résultats expérimentaux (Y)
result_effect: le groupe d’effet de l’expérience, soit, mortalité (MOR) exclusivement pour les animauxresult_conc_mean: la concentration effective moyenne en mg/Lresult_conc_mean_log: le log de la concentration effective moyenneresult_conc1_mean_binary: si la dose est moins toxique (0) ou plus toxique (1)
Variables sur les espèces (X)
tax_group: le groupe taxonomique soit,fishoucrustatax_gs: le nom de l’espèce étudiée codée commeGenre_épithètetax_lh_amd: la longévité de l’espèce en jours
Variables liées aux conditions expérimentales (X)
test_exposure_type: le type d’exposition, soit, statique (S), en flux continu d’eau (F), en flux continu d’eau avec renouvellement (R), ou non rapporté (NR)media_temperature_mean: la température moyenne du medium en degrés Celsius
Pour l’atelier, nous allons nous intéresser à modéliser les variables liées aux résultats expérimentaux. Ce seront nos variables réponse (\(Y\)). Nous testerons comment ces variables changent en fonction des variables liées aux espèces ainsi qu’aux conditions expérimentales (\(X\))
4. Exploration des données
Distribution de la concentration effective
Commençons par évaluer la distribution de la variable réponse, qui est le log de la concentration effective moyenne. Nous allons faire nos graphiques avec la librairie ggplot2.
ggplot(data = anim) +
geom_histogram(
aes(x = result_conc1_mean_log),
color = "black", fill = "#440154"
) +
ylab("Compte") +
xlab("log(concentration moyenne (mg/L))") +
theme_classic(base_size = 14) +
theme(panel.grid = element_blank())Concentration effective selon l’espèce
Vérifions si les concentrations effectives varient selon l’espèce.
ggplot(data = anim) +
geom_boxplot(
aes(
x = result_conc1_mean_log,
y = tax_gs,
fill = tax_group
),
color = "black"
) +
ylab("Espèce") +
xlab("log(concentration moyenne (mg/L))") +
scale_fill_manual(values = c("#21918c", "#440154")) +
scale_x_continuous(
breaks = seq(-8, 8, 4),
limits = c(-9, 9)
) +
theme_classic(base_size = 14) +
theme(
panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
legend.position = "top"
)Concentration effective selon l’exposition
Finalement, on peut aussi se douter que la concentration effective pourrait changer selon le type d’exposition.
ggplot(data = anim) +
geom_boxplot(
aes(
x = result_conc1_mean_log,
y = test_exposure_type,
fill = tax_group
),
color = "black"
) +
ylab("Type d'exposition") +
xlab("log(concentration moyenne (mg/L))") +
scale_fill_manual(values = c("#21918c", "#440154")) +
scale_x_continuous(
breaks = seq(-8, 8, 4),
limits = c(-9, 9)
) +
theme_classic(base_size = 14) +
theme(
panel.grid = element_blank(),
legend.position = "top"
)Sommaire de l’exploration de données
Les vérifications nous indiquent qu’il semble y avoir une structure hiérarchique dans nos données. En termes simples, cela signifie qu’il y a de la variation dans les concentrations effectives entre les différents groupes que nous avons vu, soit, les espèces ainsi que le type d’exposition.
Les modèles hiérarchiques nous permettront de modéliser cette variation afin d’en tenir compte.
5. Rappel sur les modèles linéaires
Les bases
L’objectif d’un modèle linéaire est de prédire une variable \(Y\) tout en minimisant l’écart entre les observations et les prédictions (\(\hat{Y}\)).
Au travers votre parcours, vous avez sûrement vu différents types de modèles où on utilise des variables explicatives \(X_i\) pour prédire une variable \(Y\) (ex. ANOVA, régression, ANCOVA).
Par exemple, une régression linéaire peut s’écrire avec cette équation:
\[Y = X\beta+\epsilon\]
où:
- \(Y\) est le vecteur des valeurs observées (c.-à.-d. la variable réponse)
- \(X\) est la matrice de variables explicatives (incluant l’ordonnée à l’origine)
- \(\beta\) est le vecteur de coefficients
- \(\epsilon\) est le vecteur des résidus
Visuellement, on peut représenter ce modèle avec la figure suivante:
Pour interpréter les résultats d’un modèle linéaire, il est essentiel que celui-ci respecte ces conditions d’application:
- Linéarité: la relation entre \(Y\) et \(X\) est linéaire.
- Indépendance: les observations sont indépendantes l’une de l’autre.
- Homoscédasticité: la variance des résidus est homogène le long des valeurs prédites.
- Normalité des résidus: les résidus suivent une distribution Gaussienne.
Limites des modèles linéaires
Une des limites majeures des modèles linéaires classiques est qu’ils ne peuvent s’ajuster à des variables réponse pour lesquelles les résidus suivent une distribution non-Gaussienne.
Cette situation est particulièrement commune lorsqu’on analyse des données biologiques, qui ont tendance à suivre d’autres types de distribution telle que de type Poisson ou binomiale.
Dans la section suivante, nous verrons comment les modèles linéaires généralisés nous permettent de palier à ce problème en modélisant spécifiquement la nature des résidus d’une variable \(Y\).
6. Théorie: les GLM
Dans cette section, nous verrons comment fonctionnent les modèles linéaires généralisés (GLM). Nous ferons tout d’abord un survol de trois familles de distribution. En les reconnaissant, nous pourrons modéliser nos résidus correctement lorsque nos données suivent une de ces distributions.
Pour y arriver, nous verrons aussi le concept de la fonction de lien rattachée à notre famille de distribution. La fonction de lien nous permettra d’estimer la relation linéaire entre \(Y\) et \(X\) lorsque nos résidus ne sont pas Gaussiens. Nous n’aurons donc pas à transformer nos données pour satisfaire les conditions d’application des modèles linéaires standard.
Familles de distribution
Il existe une multitude de familles de distribution qui décrivent les variables qu’on peut mesurer. Ici, l’important est de retenir que chaque distribution a des paramètres propres qui décrivent sa forme. Les GLM utilisent ces paramètres pour l’estimation des coefficients, et nous les utiliserons pour l’interprétation de nos modèles.
Voici trois exemples des familles les plus communes.
Distribution de Poisson
Certaines distributions sont discrètes et décrivent des nombres entiers exclusivement, telles que la distribution de Poisson. On voit souvent la distribution de Poisson dans les données écologiques comme l’abondance d’espèces ou le nombre d’individus d’une espèce. La distribution de Poisson est décrite par un seul paramètre, \(\lambda\), qui décrit sa moyenne et sa variance.
Voici une représentation de cette distribution avec différentes valeur de \(\lambda\):
Distribution de Bernoulli
Une autre distribution communément rencontrée dans les données écologiques est la distribution de Bernoulli. Cette distribution est décrite par le paramètre \(p\), étant la probabilité de succès (1) lorsqu’on a données de type 0 ou 1. On observe souvent cette distribution avec des données de présence-absence d’espèces, soit, typiquement lorsqu’une variable n’a que deux issues.
Cette distribution considère l’issue (0, 1) pour une seule expérience ou essai.
Voici une représentation de cette distribution avec différentes valeur de \(p\):
Distribution binomiale
Finalement, la distribution binomiale est aussi communément observée dans les données biologiques, et consiste en une extension de la distribution de Bernoulli lorsqu’on répète l’essai à plusieurs reprises.
Elle consiste donc en le nombre de succès pour un nombre fixe d’essais indépendants, où chaque essai a la même probabilité de succès \(p\). Ses paramètres sont le nombre d’essais \(n\) et la probabilité à chaque essai \(p\).
Voici une représentation de cette distribution avec différentes valeur de \(p\) et un \(n\) de 50:
Fonction de lien
Nous avons vu trois exemples de distributions souvent rencontrées avec des données biologiques.
Dans un GLM, afin que celui-ci puisse estimer la relation entre \(Y\) et \(X\), on doit donc:
- Spécifier la distribution des résidus
- Spécifier la fonction de lien pour les valeurs prédites
Cette deuxième étape consiste à linéariser la relation entre \(Y\) et \(X\). Il s’agit donc de transformer les valeurs prédites \(\hat{Y}\) pour avoir une relation linéaire.
Mathématiquement, le prédicteur linéaire dans un GLM est défini comme:
\[\eta = X\beta\]
où:
- \(\beta\) sont les coefficients
- \(X\) est la matrice des variables
La fonction de lien transforme le prédicteur linéaire à l’échelle appropriée pour \(Y\) tel que:
\[g(\mu) = \eta\]
où:
- \(\mu\) est la moyenne prédite de \(Y\) en fonction des prédicteurs \(X\)
- \(g(.)\) est la fonction de lien
Voici une table qui résume les informations relatives à la fonction de lien pour les familles que nous avons vu.
| Distribution | Nom | Fonction_de_lien | Modèle | brms |
|---|---|---|---|---|
| Gaussienne | Identité | \(g(\mu) = \mu\) | \(\mu = \mathbf{X} \boldsymbol{\beta}\) | gaussian(link="identity") |
| Poisson | Log | \(g(\mu) = \log(\mu)\) | \(\log(\mu) = \mathbf{X} \boldsymbol{\beta}\) | poisson(link="log") |
| Bernoulli | Logit | \(g(\mu) = \log\left(\frac{\mu}{1 - \mu}\right)\) | \(\log\left(\frac{\mu}{1 - \mu}\right) = \mathbf{X} \boldsymbol{\beta}\) | bernoulli(link="logit") |
| Binomial | Logit | \(g(\mu) = \log\left(\frac{\mu}{1 - \mu}\right)\) | \(\log\left(\frac{\mu}{1 - \mu}\right) = \mathbf{X} \boldsymbol{\beta}\) | binomial(link="logit") |
7. Exemple: GLM avec les données ADORE
Maintenant que nous avons vu les familles de distribution ainsi que la mécanique de la fonction de lien, nous pouvons estimer les paramètres d’un GLM avec R.
Nous allons construire un GLM simple avec une variable réponse et trois variables explicatives.
Nos variable sont:
- \(Y\): la concentration effective est 1 (plus toxique) ou 0 (moins toxique)
- \(X_1\): le type d’exposition
- \(X_2\): la température moyenne du medium degrés Celsius
- \(X_3\): la longévité de l’espèce en jours
Si on modélisait la relation entre la toxicité de la concentration effective et la température du medium avec une régression linéaire classique, on obtiendrait ceci:
Comme vous pouvez voir dans le panneau à gauche, le modèle ne décrit pas bien la relation entre le niveau de toxicité et la température. De plus, il aurait pu y avoir des cas où la droite dépasse le seuil possible de nos valeurs (ici 0 ou 1).
On peut voir à gauche que les résidus de notre modèle ne respectent pas la condition d’homoscédasticité.
Pour ces raisons, on doit donc établir une limite inférieure et supérieure à nos données en spécifiant une distribution de Bernoulli pour modéliser la toxicité de la concentration effective.
Préparation des données
Commençons par évaluer la structure des variables avant de faire le modèle.
str(anim[
, .(
result_conc1_mean_binary, test_exposure_type,
tax_lh_amd, media_temperature_mean
)
])## Classes 'data.table' and 'data.frame': 28999 obs. of 4 variables:
## $ result_conc1_mean_binary: int 0 0 0 0 0 0 0 0 0 0 ...
## $ test_exposure_type : Factor w/ 4 levels "F","NR","R","S": 3 3 3 3 3 3 3 3 3 3 ...
## $ tax_lh_amd : num 5475 5475 5475 5475 5475 ...
## $ media_temperature_mean : num 19 19 19 19 19 19 NA NA NA NA ...
## - attr(*, ".internal.selfref")=<externalptr>
Les données sont dans le format adéquat pour la modélisation.
Avans de faire le modèle, nous allons normaliser les variables explicatives continues en les ramenant à la même échelle. Ceci peut faciliter les calculs des coefficients et simplifier l’interprétation des paramètres.
# Fonction de normalisation en score Z.
standardize <- function(x) {
(x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}
# Créer deux nouvelles colones
anim[
, c("Ztax_lh_amd", "Zmedia_temperature_mean") := lapply(
.SD,
standardize
)
, .SDcols = c("tax_lh_amd", "media_temperature_mean")
]Construction du modèle
Dans cet exemple, on s’intéresse à estimer comment la toxicité de la concentration effective varie selon le type d’exposition, la longévité moyenne de l’espèce, et la température du médium durant l’expérience.
Voici l’équation de notre modèle:
\[\eta = X\beta + Z\gamma + 1\beta_0 + \epsilon\]
où:
- \(\eta\) est le prédicteur linéaire
- \(X\) est la matrice de conception du type d’exposition
- \(\beta\) est le vecteur des coefficients du type d’exposition
- \(Z\) est la matrice de conception des prédicteurs linéaires (longévité et température)
- \(\gamma\) est le vecteur des coefficients du type d’exposition des prédicteurs linéaires
- \(1\) est le vecteur de 1 pour l’ordonnée à l’origine
- \(\beta_0\) est l’ordonnée à l’origine pour le niveau de référence
La fonction de lien du modèle est:
\[logit(p) = log(\frac{p}{1-p})\]
où:
\(p\) est la probabilité que la concentration effective soit toxique (1) pour chaque observation.
On estime \(p\) en utilisant l’inverse de la fonction logit:
\[p = \text{logit}^{-1}(\eta) = \frac{1}{1 + e^{-\eta}}\]
Spécifions ce modèle dans brms avec la fonction brmsformula():
form1 <- brmsformula(
result_conc1_mean_binary ~
test_exposure_type +
Ztax_lh_amd +
Zmedia_temperature_mean
)Informations à priori
Comme nous l’avons vu dans la section sur l’inférence Bayésienne, la distribution postérieure est estimée en multipliant la distribution à priori et la distribution de vraisemblance des données (les preuves).
Par défaut, brms emploie des distributions à priori très larges qui ne réflètent pas nos hypothèses/connaissances sur notre système. Définir les connaissances à priori est donc très important en inférence Bayésienne, et je vous suggère fortement de faire ce processus. C’est ce qui différencie l’inférence Bayésienne de l’approche fréquentiste.
Ici, puisque je ne connais pas assez bien l’écotoxicologie, je spécifie des “priors” relativement informatifs suivant une distribution Gaussienne avec une moyenne de 0 et un écart-type de 1.
Cela signifie que je crois que les valeurs des paramètres devraient se situer autour de 0 avec des valeurs positives allant de +1 à -1 écart-types.
Les informations à priori à spécifier sont:
- sur l’ordonnée à l’origine globale (la moyenne du type d’exposition
F) - sur les différences de moyenne par rapport à
F - sur les pentes pour la température et la longévité
Pour vous aider à comprendre, imaginez que vous savez avec la littérature que la dose effective devrait être plus toxique avec une augmentation de la température. Ainsi, pour la pente de ce paramètre, vous pourriez spécifier une distribution Gaussienne avec une moyenne de 1 et un écart-type de 0.5. Vous diriez donc à votre modèle que vous assumez que cette valeur devrait être positive, se situant autour de 1, avec peu de variation autour de cette valeur.
priors1 <- c(
# Informations à priori sur l'ordonnée à l'origine
set_prior(
"normal(0, 1)",
class = "Intercept"
),
# Informations à priori sur les différences de moyenne
set_prior(
"normal(0, 1)",
class = "b",
coef = c(
"test_exposure_typeNR",
"test_exposure_typeR",
"test_exposure_typeS"
)
),
# Informations à priori sur les pentes
set_prior(
"normal(0, 1)",
class = "b",
coef = c(
"Zmedia_temperature_mean",
"Ztax_lh_amd"
)
)
)Estimation des paramètres
Spécifications pour l’estimation MCMC
À titre de rappel, brms utilise l’algorithme de HMC-NUTS pour estimer les paramètres du modèle.
L’agorithme va effectuer un nombre \(n\) d’itérations que nous allons spécifier. À chaque itération, l’algorithme estimera la valeur des paramètres d’intérêt de notre modèle.
Généralement, on effectue \(n\) itérations sur 4 chaînes de Markov qu’on combine pour avoir les distributions postérieures. Pour une bonne précision, on voit souvent des distributions postérieures avec \(\ge 1000\) échantillons au total.
Voici les spécifications pour l’estimation des paramètres par HMC-NUTS. Pour chacune des 4 chaînes, nous aurons:
- n: 2000
- réchauffement: n / 2
- itérations: n - réchauffement
- échantillonnage: itérations / échantillonnage = 1000
Notre modèle fera donc 1000 réchauffements (warmup = 1000), 1000 itérations (iter = 2000 où 2000-1000), et échantillonera une valeur de paramètre à chaque itération (thin = 1).
Nous obtiendrons ainsi 1000 échantillons par chaîne que nous combinerons pour un total de 4000 échantillons postérieurs par paramètre.
Accélérer l’estimation (facultatif)
Il est à noter qu’estimer un modèle par MCMC peut prendre beaucoup de temps, particulièrement lorsqu’il y a beaucoup de données et qu’on veut faire beaucoup d’itérations.
Pour cette raison, vous pouvez tirer profit des coeurs du processeur de votre ordinateur pour accélérer l’estimation.
Nous allons donc vérifier le nombre de coeurs que nous avons.
parallel::detectCores()## [1] 16
La fonction indique que j’ai 16 coeurs disponibles dans mon ordinateur. On doit toujours en réserver un pour les processus internes.
Si vous en avez seulement 4 ou moins, passez à la section “Si vous avez \(\le4\) coeurs”
Si vous avez \(\ge 6\) coeurs
Je vais utiliser 12 coeurs qui effectueront un processus parallèle sur 4 chaînes.
Je spécifie 12 et non 14 car 14 / 4 = 3.5, un chiffre impair. On doit toujours avoir un nombre de coeurs pairs pour effectuer les tâches en parallèle.
Puisque j’utiliserai 12 coeurs en parallèle sur 4 chaînes, 3 coeurs seront sollicités par chaîne. Ceci est spécificé dans l’argument threads = threading(3).
Si vous avez \(\le 4\) coeurs
Si vous avez 4 coeurs ou moins sur votre ordinateur, vous ne pourrez pas effectuer le modèle en parallèle. Vous devrez ainsi estimer les paramètres un échantillon et une chaîne à la fois, ce qui prendra plus de temps.
Toutefois, il est possible que 3 chaînes effectuent les calculs simultanément si vous avez 4 coeurs. Vous pouvez ainsi remplacer l’argument threads = threading(3) par cores = parallel::detectCores()-1. Vous utiliserez donc 3 coeurs pour que 3 chaînes puissent échantilloner les paramètres simultanément, un échantillon à la fois.
Ultimement, si vous avez un gros modèle à faire rouler, je suggère de l’effectuer sur un serveur de calcul si vous y avez accès.
Exécution de l’estimateur
Nous sommes donc prêt.e.s à estimer les paramètres:
glm1 <- brm(
# La formule spécifiée
formula = form1,
# La famille de distribution
family = bernoulli(link = "logit"),
# Les paramètres d'itérations pour le MCMC
iter = 2000,
warmup = 1000,
thin = 1,
chains = 4,
# Paramètres pour accélérer les calculs
# Utiliser seulement si vous avez >4 coeurs
threads = threading(3),
# Si vous avez 4 coeurs ou moins
# cores = parallel::detectCores()-1,
# Estimateur
backend = "cmdstanr",
# Reproductibilité
seed = 123,
# Paramètres de contrôle
control = list(adapt_delta = 0.95),
# Distributions à priori
prior = priors1,
# Échantillonage des distributions à priori
# Ceci peut être utile à des fins de comparaison
sample_prior = TRUE,
# Données
data = anim
)Puisque l’estimation peut prendre quelques temps, on créé un dossier où nous allons sauvegarder le modèle.
chemin <- file.path(getwd(), "outputs")
dir.create(chemin)On sauvegarde le modèle en format .rds qui pourra facilement être réimporté dans la session au besoin.
saveRDS(
object = glm1,
file = file.path(chemin, "glm1.rds")
)Vérifications du modèle
Maintenant que notre modèle a terminé son exécution, il est essentiel de faire des inspections pour s’assurer que les paramètres estimés sont fiables. Ces étapes d’inspections sont cruciales pour s’assurer que nos interprétations biologiques font du sens. Autrement, nous pourrions avoir de mauvaises surprises.
Commençons par extraire les échantillons des distributions à posteriori dans un data.frame pour les manipuler.
params <- as_draws_df(
x = glm1, add_chain = TRUE,
regex = TRUE
)Convergence des chaînes
La première chose que nous allons faire est de vérifier si les chaînes de Markov ont bien convergé pour chacun de nos paramètres d’intérêt. Nous allons faire ceci en inspectant des “trace plots”.
Nous pouvons déduire que les chaînes ont convergé si les “trace plots” sont homogènes.
bayesplot::mcmc_trace(
params,
regex_pars = "b",
np = nuts_params(glm1)
)Les chaînes de notre modèle semblent avoir très bien convergé.
Une autre étape d’inspection de convergence des chaînes de Markov est d’évaluer les \(\hat{R}\) et les tailles d’échantillon effectives pour chaque paramètre d’intérêt.
Dans Bürkner (2017), on peut lire que la taille d’échantillon effective est: “le nombre d’échantillons indépendants de la distribution postérieure qui serait attendu pour produire la même erreur standard de la moyenne postérieure que celle obtenue à partir des échantillons dépendants retournés par l’algorithme MCMC”. Dans
brms, le résumé du modèle fournit une estimation des tailles d’échantillon effectives de la partie centrale et des queues de la distribution postérieure. LeBulk_ESSrapporte l’efficacité de l’échantillonnage de la partie centrale de la distribution postérieure (c’est-à-dire la médiane ou la moyenne), tandis que leTail_ESSdiagnostique l’efficacité de l’échantillonnage des queues de la distribution postérieure (c’est-à-dire les quantiles 5% et 95%). Nous utilisons un seuil de <100 comme règle de décision pour savoir si les chaînes ont convergé pour les paramètres (Vehtari et al. 2021).En revanche, le \(\hat{R}\) est un outil de diagnostic quantitatif courant pour savoir si les chaînes ont convergé (Gelman et al. 2013; Bürkner 2017). Il compare la variance à l’intérieur d’une chaîne à la variance entre les chaînes. Différentes versions et seuils ont été proposés pour cette mesure (Gelman and Rubin 1992; Gelman et al. 2013; Vehtari et al. 2021), mais nous utiliserons la plus parcimonieuse suggérée par Vehtari et al. (2021) avec un seuil fixé à <1.01.
Dans brms, le \(\hat{R}\), le Bulk_ESS, et le Tail_ESS apparaissent toujours à côté des moyennes postérieures des paramètres affichées avec la fonction summary(). Examinons les valeurs pour nos paramètres:
glm1_t <- data.table(
round(summary(glm1, prob = 0.89, robust = TRUE)$fixed, digits = 3),
keep.rownames = TRUE
)
setnames(glm1_t, old = "rn", new = "Parameter")
glm1_t[, c(1, 6:8)]Les valeurs sont très bien!
Ajustement du modèle
Nous allons maintenant examiner si notre modèle a bien prédit les données en utilisant les vérifications prédictives postérieures. Les vérifications prédictives postérieures utilisent des données simulées à partir de la distribution prédictive postérieure d’un modèle et les comparent aux données brutes utilisées pour ajuster le modèle. Gabry et al. (2019) suggèrent qu’un modèle bien ajusté devrait générer des données ressemblant aux données qui ont été observées (bruts). Ici, \(y\) est la distribution de nos données observées, et \(y_{rep}\) est la distribution des données simulées à partir de la distribution prédictive postérieure du modèle.
La production des graphiques peut prendre du temps en fonction de la complexité du modèle et de la taille des données. Vous pouvez spécifier le nombre d’échantillons à afficher en utilisant l’argument ndraws dans la fonction pp_check() pour accélérer le processus.
pp_check(glm1, ndraws = 100)Notre modèle semble très bien ajusté aux données!
En somme, nos inspections indiquent que le modèle est très bon. Toutefois, je vous encourage à consulter la section sur les inspections additionnelles pour mieux diagnostiquer vos modèles.
Interprétations
Sommaire des paramètres
Par défaut, la fonction summary() dans brms rapporte la moyenne de la distribution postérieure avec les intervales de compatibilité inférieurs et supérieurs à 95%.
Toutefois, il est suggéré en analyse Bayésienne de rapporter les intervalles à 89% car ils procurent davantage de stabilité. De plus, la médiane est une meilleure mesure de tendance centrale lorsque la distribution postérieure n’est pas symétrique.
Nous avons donc ajouté des arguments à la fonction summary précédemment en faisant ce code:
glm1_t <- data.table(
round(summary(glm1, prob = 0.89, robust = TRUE)$fixed, digits = 3),
keep.rownames = TRUE
)
setnames(glm1_t, old = "rn", new = "Parameter")où robust = TRUE indique de rapporter la médiane, et prop=0.89 les intervalles à 89%. Au final, la décision de l’intervalle et arbitraire, et le mieux est de visualiser la distribution entière de l’estimé. Nous ferons ceci plus tard.
Visualisons la table de sommaire:
knitr::kable(glm1_t[, c(1, 2, 4, 5)])| Parameter | Estimate | l-89% CI | u-89% CI |
|---|---|---|---|
| Intercept | -0.404 | -0.459 | -0.348 |
| test_exposure_typeNR | 0.449 | 0.353 | 0.547 |
| test_exposure_typeR | -0.027 | -0.127 | 0.070 |
| test_exposure_typeS | 0.161 | 0.097 | 0.224 |
| Ztax_lh_amd | 0.077 | 0.055 | 0.099 |
| Zmedia_temperature_mean | -0.133 | -0.155 | -0.110 |
On peut voir que:
- Il y a des différences de toxicité selon le type d’exposition
- La probabilité que la CL50 soit plus toxique augmente avec la longévité
- La probabilité que la CL50 soit plus toxique diminue avec la température moyenne
Nous allons visualiser ces résultats avec des figures.
Figures
Calcul des prédictions
Tout d’abord, on extrait les prédictions pour la longévité moyenne ainsi que la température moyenne avec la fonction conditional_effects() qui est spécifique à brms.
Quelques spécifications:
- l’argument
spaghetti=va chercher des prédictions posterieures - l’argument
ndraws=spécifie le nombre de prédictions postérieures (ici, 100) - pour chacune des 100 droites, nous aurons 100 valeurs. C’est le paramètre par défaut de
conditional_effects()
Il est possible de tracer le graphique directement à partir de la fonction. Toutefois, nous allons plutôt conserver les valeurs dans une table afin d’avoir plus de flexibilité sur les paramètres graphiques.
Température
tab1 <- conditional_effects(
x = glm1, method = "posterior_epred",
effects = "Zmedia_temperature_mean",
spaghetti = TRUE, ndraws = 100
)
spag1 <- attr(tab1$Zmedia_temperature_mean, "spaghetti")Longévité
tab2 <- conditional_effects(
x = glm1, method = "posterior_epred",
effects = "Ztax_lh_amd",
spaghetti = TRUE, ndraws = 100
)
spag2 <- attr(tab2$Ztax_lh_amd, "spaghetti")Type d’exposition
Ici, on va utiliser les distributions postérieures des moyennes prédites de chaque type d’exposition. Nous avons déjà extrait ceci dans une étape ultérieure et assigné à l’objet params.
Cela nous permettra de comparer les moyennes entre chaque type d’exposition avec une approche Bayésienne.
Visualisons l’objet:
paramsNous allons utiliser les 4 premières colones pour produire le graphique. Ces 4 colones sont:
b_Intercept: la moyenne de l’exposition de typeFb_test_exposure_typeNRest la différence entreNRetFb_test_exposure_typeRest la différentre entreRetFb_test_exposure_typeSest la différentre entreSetF
Nous devons manipuler un peu la table pour faire le graphique.
On veut:
- calculer la moyenne pour les expositions
NR,R, etS - ramener les valeurs à une échelle de probabilité
- calculer la médiane des distributions postérieures
Commençons par l’étape 1.
# Transformer en data.table pour manipuler
paramsdt <- data.table(params)
# On calcule les vraies moyennes
paramsdt[, c(2:4) := paramsdt[, c(2:4)] + b_Intercept]Ensuite l’étape 2.
Ici, souvenez-vous que la fonction de lien nous a permis d’estimer un modèle sur une échelle linéaire. Toutefois, cette échelle est difficilement interprétable. On doit donc ramener les prédictions à l’échelle de notre variable, soit, en probabilités.
# Changer en format long
preds3 <- melt(
data = paramsdt[, c(1:4)],
variable.name = "exposure_type",
measure = colnames(paramsdt[, c(1:4)])
)
# Ramener à l'échelle de probabilité
preds3[, value := plogis(value)]Finalement, l’étape 3.
# Calculer les médianes
medians <- preds3[, .(median = median(value)), by = exposure_type]Visualisation
On peut maintenant produire les graphiques. À titre de rappel, on a la même pente pour l’ensemble des types d’exposition, et des ordonnées à l’origine qui diffèrent (c-à-d, les moyennes de toxicité diffèrent). Ici, la droite correspond au type d’exposition F.
Commençons avec les différences de moyennes.
ggplot(
data = preds3,
aes(x = value, fill = exposure_type)
) +
geom_density(alpha = 0.5, colour = "black") +
geom_vline(
data = medians,
aes(xintercept = median, color = exposure_type),
linetype = "dashed"
) +
scale_fill_manual(
labels = c(
"Flux", "NS",
"Renouvellement", "Statique"
),
values = c("#440154", "#3d4d8a", "#24878e", "#40bd72")
) +
scale_colour_manual(
labels = c(
"Flux", "NS",
"Renouvellement", "Statique"
),
values = c("#440154", "#3d4d8a", "#24878e", "#40bd72")
) +
#scale_x_continuous(breaks = seq(0.25, 1, 0.25), limits = c(0.25, 1)) +
xlab("\nProbabilité d'être plus toxique") +
ylab("Densité\n") +
labs(fill = "Exposition", colour = "Exposition") +
theme_classic(base_size = 14) +
theme(legend.position = "top")Finalement, les droites pour la température et la longévité pour le type d’exposition F.
preds1 <- ggplot() +
geom_point(
data = anim,
aes(
y = result_conc1_mean_binary,
x = Zmedia_temperature_mean
),
alpha = 0.05
) +
geom_line(
data = spag1,
aes(
x = Zmedia_temperature_mean,
y = estimate__,
group = sample__
),
linewidth = 1, alpha = 0.05,
colour = "#440154"
) +
ylab("Toxicité\n") +
scale_y_continuous(
breaks = seq(0, 1, 0.25),
limits = c(0, 1)
) +
xlab("\nTempérature moyenne normalisée") +
theme_classic(base_size = 14)
preds2 <- ggplot() +
geom_point(
data = anim,
aes(
y = result_conc1_mean_binary,
x = Ztax_lh_amd
),
alpha = 0.05
) +
geom_line(
data = spag2,
aes(
x = Ztax_lh_amd,
y = estimate__,
group = sample__
),
linewidth = 1, alpha = 0.05,
colour = "#440154"
) +
ylab("Toxicité\n") +
scale_y_continuous(
breaks = seq(0, 1, 0.25),
limits = c(0, 1)
) +
xlab("\nLongévité moyenne") +
theme_classic(base_size = 14)
ggarrange(preds1, preds2)8. Théorie: les HGLM
Complexité des données biologiques
En explorant nos données, nous avons vu qu’il semble y avoir une structure hiérarchique, où la concentration effective varie entre les espèces.
Ceci est un exemple concret démontrant comment les données biologique sont souvent complexes, et structurées à différentes échelles.
Reprenons l’exemple de la concentration effective en fonction de la température du médium. Cette fois, on ajoute les espèces dans la figure. Ici, chaque couleur correspond à un taxon différent.
À partir de ce graphique, on peut se poser deux questions:
- Est-ce que les taxons ont la même concentration effective moyenne?
- Est-ce que la relation entre la concentration effective et la température change selon les taxons?
Pour la première question, notre réflexe pourrait être de faire le même GLM que nous avons fait à la section précédente, et d’ajouter le taxon comme variable catégorique.
Pour la seconde question, on pourrait spécifier un modèle de type ANCOVA, en modélisant une interaction entre la température et le taxon.
Pour répondre aux deux questions, on aurait ce modèle:
result_conc1_mean_log ~
test_exposure_type +
Ztax_lh_amd +
Zmedia_temperature_mean +
tax_gs +
Zmedia_temperature_mean:tax_gsToutefois, on se retrouve avec un gros problème ici. Visualisons-le:
Nous avons 105 espèces, et donc 105 droites à estimer. Cela signifie que le modèle doit estimer:
- 105 ordonnées à l’origine
- 105 pentes
C’est donc 210 paramètres de + à estimer! Si on utilisait l’approche fréquentiste pour l’analyse de données, ça nous coûterait cher en degrés de liberté. Ceci est moins problématique pour l’approche Bayésienne car on n’utilise pas la valeur p.
Fréquentiste ou Bayésien, le plus important, c’est qu’en estimant 210 paramètres, la sortie de la fonction summary() serait cauchemardesque. Pensez-y! Plus haut, nous avions un Intercept_b qui utilisait le type d’exposition F comme référence. Ensuite, le modèle calculait la différence entre les autres types d’exposition et F. Je vous laisse imaginer celà avec 105 taxons…
L’avantage des modèles hiérarchiques
Les modèles hiérarchiques nous permettent de palier aux problèmes que nous avons vu plus haut:
- permettent de séparer des effets populationnels d’effets groupés
- permettent de modéliser la structure hiérarchique de nos données en réduisant le nombre de paramètres à estimer
- exploitent les paramètres populationnels pour estimer les paramètres de groupe même quand la taille d’échantillon est faible pour un groupe
De plus, rappelez vous qu’une des conditions de base dans les modèles linéaires est que les échantillons doivent être indépendants. Dans l’exemple que nous avons vu, nous avons plusieurs observations par espèce, ce qui contrevient à cette condition d’indépendance.
Les modèles hiérarchiques permettent donc de modéliser des facteurs associés à la structure de vos expériences. Voici quelques exemples:
- On peut intégrer les quadrats, ou les sites
- On peut ajouter les individus si des données sont prises à répétition
- Permettent de considérer l’aspect temporel dans les données
Nous verrons comment ces principes fonctionnent dans la prochaine section avec les concepts d’effets fixes et d’effets aléatoires.
Effets fixes et aléatoires
Reprenons l’exemple du modèle que nous avons construit précédemment.
On s’intéresse toujours à estimer comment la concentration effective varie selon le type d’exposition, la longévité moyenne de l’espèce, et la température du médium durant l’expérience. On ajoute toutefois les espèces comme effet aléatoire.
Ici, on veut donc savoir si les espèces varient dans leur concentration effective moyenne, et si la relation entre celle-ci et la température change d’une espèce à l’autre.
Voici l’équation de notre modèle:
\[\eta = X\beta + Z\gamma + U_1\delta_1 + U_2\delta_2 + 1\beta_0 + \epsilon\]
où les paramètres ajoutés sont:
- \(U_1\) est la matrice de conception pour les ordonnées à l’origine aléatoires
- \(\delta_1\) est le vecteur des ordonnées à l’origine aléatoires des espèces
- \(U_2\) est la matrice de conception pour les pentes aléatoires
- \(\delta_2\) est le vecteur des pentes aléatoires des espèces
Les effets fixes, représentés par les coefficients \(\beta\) et \(\gamma\), modélisent les effets moyens des variables explicatives sur la variable réponse. Ils sont constants pour toutes les unités observées (c.-à.-d., les espèces par exemple). Ce sont les effets que nous avons estimés précédemment dans notre GLM. Celà ne change donc pas. Les effets fixes peuvent être des variables catégoriques ou des variables continues
Les effets aléatoires, quant à eux, représentés par les coefficients \(\delta_1\) et \(\delta_2\), modélisent la variation entre les différents niveaux d’une variable catégorique. Dans le contexte de notre exemple, on intègre les différences entre les espèces en termes de sensibilité à la concentration en ajoutant un effet aléatoire. Les effets aléatoires sont strictement des variables catégoriques
En ajoutant cet effet aléatoire, le modèle peut capturer une part supplémentaire de la variation de \(Y\) qui ne peut pas être expliquée par les effets fixes. Cela permet de mieux ajuster le modèle aux données et d’améliorer les prédictions.
Les types de modèles hiérarchiques
Il existe deux types de modèles hiérarchiques correspondant aux effets aléatoires estimés.
Pour le premier type, on assume que la relation entre \(Y\) et \(X\) est la même pour l’ensemble de nos groupes, mais que la moyenne de \(Y\) varie entre les groupes. C’est donc un modèle où l’on estime une moyenne pour chacun des niveaux de notre facteur d’intérêt, mais pour lequel la relation entre \(Y\) et \(X\) est la même que l’effet fixe pour l’ensemble des niveaux du facteur. Ce type de modèle hiérarchique se nomme modèle d’ordonnées à l’origine aléatoires.
Le second type de modèle est construit à partir de la mécanique du modèle d’ordonnées à l’origine aléatoires. C’est le modèle de pentes aléatoires. Comme le nom l’indique, c’est un modèle où chaque niveau du facteur peut maintenant avoir sa propre relation entre \(Y\) et \(X\).
Comme nous l’avons soulevé plus haut, un avantage majeur des modèles hiérarchiques est qu’ils réduisent le nombre de paramètres à estimer. Au début de cette section, nous avons vu que si nous faisions une ANCOVA, nous aurions 210 paramètres à estimer. En exploitant les effets aléatoires nous réduisons ce nombre de paramètres à 2 seulement!
La raison est qu’au lieu d’estimer la différence d’ordonnée à l’origine pour chaque niveau du facteur ainsi que la différence de pente, un modèle hiérarchique va plutôt calculer chaque valeur pour chaque niveau (c.-à-d., moyenne et pente), et ensuite estimer une variance.
On a donc:
- la variance entre les ordonnées à l’origine (c.-à-d., les moyennes des niveaux du facteur)
- la variance entre les pentes (c.-à-d., la relation entre \(Y\) et \(X\) pour chaque niveau du facteur)
Ce que la variance nous dit:
- une variance faible entre les ordonnées à l’origine = la moyenne de \(Y\) de chaque espèce est semblable
- une variance faible entre les pentes = la relation entre \(Y\) et \(X\) pour chaque espèce est semblable
On assume que les ordonnées à l’origine et les pentes suivent une distribution normale telle que:
\[ \begin{bmatrix} \delta_1 \\ \delta_2 \end{bmatrix} \sim \text{MVN} \left( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \Omega \right) \]
où \(\Omega\) est une matrice de covariance telle que:
\[ \Omega = \begin{bmatrix} \text{var}(\delta_1) & \text{cov}(\delta_1, \delta_2) \\ \text{cov}(\delta_2, \delta_1) & \text{var}(\delta_2) \end{bmatrix} \]
Ordonnées à l’origine aléatoires
Voici à quoi ressemble ce type de modèle pour des données hiérarchiques.
On peut voir que :
- chaque groupe a une ordonnée à l’origine différente
- tous les groupes on la même pente, suivant la tendance populationnelle (droite pointillée)
- pour l’effet aléatoire, nous avons 1 paramètre étant la variance entre les ordonnées à l’origine
Pentes aléatoires
Voici à quoi ressemble ce type de modèle pour des données hiérarchiques.
On peut voir que :
- chaque groupe a une ordonnée à l’origine différente
- chaque groupe a une pente différente
- pour l’effet aléatoire, nous avons 2 paramètres étant la variance entre les ordonnées à l’origine et la variance entre les pentes
9. Exemple: HGLM avec les données ADORE
Pour mieux visualiser les relations, nous ferons cette fois-ci un HLM en utilisant le log de la concentration effective moyenne.
form2 <- brmsformula(
result_conc1_mean_log ~
test_exposure_type +
Ztax_lh_amd +
Zmedia_temperature_mean +
(1 + Zmedia_temperature_mean | tax_gs)
)priors2 <- c(
# Informations à priori sur l'ordonnée à l'origine
set_prior(
"normal(0, 1)",
class = "Intercept"
),
# Informations à priori sur les différences de moyenne
set_prior(
"normal(0, 1)",
class = "b",
coef = c(
"test_exposure_typeNR",
"test_exposure_typeR",
"test_exposure_typeS"
)
),
# Informations à priori sur les pentes
set_prior(
"normal(0, 1)",
class = "b",
coef = c(
"Zmedia_temperature_mean",
"Ztax_lh_amd"
)
),
# Informations à priori sur les ordonnées à l'origine aléatoires
set_prior(
"normal(0, 1)",
class = "sd",
coef = "Intercept",
group = "tax_gs"
),
# Informations à priori sur les pentes aléatoires
set_prior(
"normal(0, 1)",
class = "sd",
coef = "Zmedia_temperature_mean",
group = "tax_gs"
)
)hglm <- brm(
# La formule spécifiée
formula = form2,
# La famille de distribution
family = gaussian(),
# Les paramètres d'itérations pour le MCMC
iter = 2000,
warmup = 1000,
thin = 1,
chains = 4,
# Paramètres pour accélérer les calculs
# Utiliser seulement si vous avez >4 coeurs
threads = threading(3),
# Si vous avez 4 coeurs ou moins
# cores = parallel::detectCores()-1,
# Estimateur
backend = "cmdstanr",
# Reproductibilité
seed = 123,
# Paramètres de contrôle
control = list(adapt_delta = 0.95),
# Distributions à priori
prior = priors2,
# Échantillonage des distributions à priori
# Ceci peut être utile à des fins de comparaison
sample_prior = TRUE,
# Données
data = anim
)On sauvegarde le modèle en format .rds qui pourra facilement être réimporté dans la session au besoin.
saveRDS(
object = hglm,
file = file.path(chemin, "hglm.rds")
)Informations complémentaires
Inspections additionnelles
While the model verifications that we just did are probably among the most important ones, they cover only a small portion of the checks that you should perform to ensure that your model has been correctly specified. To go further, here is a list of useful links for additional checks:
- MCMC diagnostics with the
bayesplotpackage - Additional posterior predictive checks with the
bayesplotpackage - Additional residuals inspections with the
tidybayespackage. For instance, you could inspect the relationship between the residuals and your fixed effects to see if the model failed at describing some features. This will be important for biological interpretations. - Leave-one-out cross-validation using Pareto-smoothed importance sampling for predictive performance and model comparisons. For important references consult the loo FAQ, the roaches model example, and the loo package vignette
Moreover, note that MCMC estimation can come with different convergence problems. Those problems can arise from different processes, but brms (and other packages using MCMC) should print warnings outlining the potential problems after the model has finished running. Here is a useful link describing each of those problems and their solutions.
References
Aczel, Balazs, Rink Hoekstra, Andrew Gelman, Eric-Jan Wagenmakers, Irene G. Klugkist, Jeffrey N. Rouder, Joachim Vandekerckhove, et al. 2020. “Discussion Points for Bayesian Inference.” Nature Human Behaviour 4 (6): 561–63. https://doi.org/10.1038/s41562-019-0807-z.
Bürkner, Paul-Christian. 2018. “Advanced Bayesian Multilevel Modeling with the R Package Brms.” The R Journal 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.
———. 2017. “Brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.
Gabry, Jonah, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman. 2019. “Visualization in Bayesian Workflow.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 182 (2): 389–402. https://doi.org/10.1111/rssa.12378.
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. Third. New York: Chapman and Hall/CRC. https://doi.org/10.1201/b16018.
Gelman, Andrew, and Donald B. Rubin. 1992. “Inference from Iterative Simulation Using Multiple Sequences.” Statistical Science 7 (4). https://doi.org/10.1214/ss/1177011136.
Kruschke, John K. 2021. “Bayesian Analysis Reporting Guidelines.” Nature Human Behaviour 5 (10): 1282–91. https://doi.org/10.1038/s41562-021-01177-7.
McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and STAN. Second. New York: Chapman and Hall/CRC. https://doi.org/10.1201/9780429029608.
Piironen, Juho, and Aki Vehtari. 2017. “Comparison of Bayesian Predictive Methods for Model Selection.” Statistics and Computing 27 (3): 711–35. https://doi.org/10.1007/s11222-016-9649-y.
Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and WAIC.” Statistics and Computing 27 (5): 1413–32. https://doi.org/10.1007/s11222-016-9696-4.
Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. “Rank-Normalization, Folding, and Localization: An Improved $\widehat{}R{}$ for Assessing Convergence of MCMC (with Discussion).” Bayesian Analysis 16 (2): 667–718. https://doi.org/10.1214/20-BA1221.
Reproductibilité
Voici les informations minimales de mon système pour reproduire les analyses avec les mêmes versions de R ainsi que les librairies.
sessionInfo()## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Pop!_OS 22.04 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_CA.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_CA.UTF-8 LC_COLLATE=en_CA.UTF-8
## [5] LC_MONETARY=en_CA.UTF-8 LC_MESSAGES=en_CA.UTF-8
## [7] LC_PAPER=en_CA.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_CA.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lme4_1.1-34 Matrix_1.6-1.1 knitr_1.46 viridis_0.6.4
## [5] viridisLite_0.4.2 ggpubr_0.6.0 ggplot2_3.4.4 bayesplot_1.10.0
## [9] brms_2.20.4 Rcpp_1.0.12 data.table_1.14.8
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.6 colorspace_2.1-0 ggsignif_0.6.4
## [4] estimability_1.4.1 markdown_1.11 QuickJSR_1.0.7
## [7] base64enc_0.1-3 farver_2.1.1 rstan_2.32.3
## [10] DT_0.30 fansi_1.0.6 mvtnorm_1.2-3
## [13] splines_4.1.2 bridgesampling_1.1-2 codetools_0.2-18
## [16] cachem_1.1.0 shinythemes_1.2.0 jsonlite_1.8.8
## [19] nloptr_2.0.3 broom_1.0.6 shiny_1.8.1.1
## [22] compiler_4.1.2 emmeans_1.8.8 backports_1.4.1
## [25] fastmap_1.2.0 cli_3.6.2 later_1.3.2
## [28] htmltools_0.5.8.1 prettyunits_1.2.0 tools_4.1.2
## [31] igraph_1.5.1 coda_0.19-4.1 gtable_0.3.4
## [34] glue_1.7.0 reshape2_1.4.4 dplyr_1.1.4
## [37] posterior_1.5.0 carData_3.0-5 jquerylib_0.1.4
## [40] vctrs_0.6.5 nlme_3.1-155 crosstalk_1.2.0
## [43] tensorA_0.36.2 xfun_0.44 stringr_1.5.1
## [46] ps_1.7.5 mime_0.12 miniUI_0.1.1.1
## [49] lifecycle_1.0.4 gtools_3.9.4 rstatix_0.7.2
## [52] MASS_7.3-55 zoo_1.8-12 scales_1.2.1
## [55] colourpicker_1.3.0 promises_1.3.0 Brobdingnag_1.2-9
## [58] parallel_4.1.2 inline_0.3.19 shinystan_2.6.0
## [61] yaml_2.3.8 curl_5.1.0 gridExtra_2.3
## [64] loo_2.6.0 StanHeaders_2.26.28 sass_0.4.9
## [67] stringi_1.8.4 highr_0.10 dygraphs_1.1.1.6
## [70] checkmate_2.3.0 boot_1.3-28 pkgbuild_1.4.2
## [73] cmdstanr_0.6.1 rlang_1.1.3 pkgconfig_2.0.3
## [76] matrixStats_1.0.0 distributional_0.3.2 evaluate_0.23
## [79] lattice_0.20-45 purrr_1.0.2 rstantools_2.3.1.1
## [82] htmlwidgets_1.6.2 labeling_0.4.3 cowplot_1.1.1
## [85] processx_3.8.2 tidyselect_1.2.1 plyr_1.8.9
## [88] magrittr_2.0.3 bookdown_0.37 R6_2.5.1
## [91] generics_0.1.3 pillar_1.9.0 withr_3.0.0
## [94] mgcv_1.9-0 xts_0.13.1 abind_1.4-5
## [97] tibble_3.2.1 crayon_1.5.2 car_3.1-2
## [100] utf8_1.2.4 rmarkdown_2.25 grid_4.1.2
## [103] callr_3.7.3 rmdformats_1.0.4 threejs_0.3.3
## [106] digest_0.6.35 xtable_1.8-4 tidyr_1.3.1
## [109] httpuv_1.6.15 RcppParallel_5.1.7 stats4_4.1.2
## [112] munsell_0.5.0 bslib_0.7.0 shinyjs_2.1.0